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We present the results of an experimental and numerical investigation of a turbulent flow over 
a backward-facing step in a channel. Experimental data are visualized using a Particle Image 
Velocimetry (PIV) device. As a mathematical model we used quasi-hydrodynamic (QHD) equa- 
tions. We have carried out numerical modeling of the flow for three configurations that were also 
studied experimentally The computed flow proved to be non-stationary. The averaged flow in the 
separation zone is in a good agreement with the experimental data. 

1 Introduction 

In this paper we compare the results of the experimental and numerical investigation of a turbulent 
flow with medium Reynolds numbers over a backward-facing step in a channel. 

An adequate description of a non-stationary flow in the separation zone behind the backward 
step is a considerable problem both in experimental and numerical aspects. As far as the numerical 
methods are concerned, a great variety of approaches are currently proposed to describe a turbu- 
lent flow. From the experimental point of view, visualization is a difficult task due to the great 
differences in the velocity of the flow in the main stream and in the separation zone, which ham- 
pers simultaneous visualization of the flow above the step and in the separation zone. Additional 
dificulties arise from the nonstationary character of the flow. 

At present there is a significant number of computations of flows over a backward-facing step in 
a channel using the averaged Navier-Stokes equations (see pQ - jS] and literature cited there). The 
length of the separation zone related to the height of the step L s /h in turbulent regime practically 
does not depend on the Reynolds number and varies - according to estimations pQ - [S] - from 5 
to 8 depending on the geometry of the problem. 

Averaged equations aimed at obtaining stationary numerical solutions, but the ambiguity in 
choosing the parameters of the model may lead to non-stationary flows behind a step as well [H]- 
Meanwhile even in the case of a developed turbulent flow the question of validity of the turbulence 
model remains. Non-stationary turbulent flow regimes may be obtained using the LES (Large- 
Eddy Simulation) or DNS (Direct Numerical Simulation) approaches, e.g. jH] - [H]. At present 
there is a number of results in numerical 3D-simulation of flows over a backward-facing step (e.g. 

The experiments described in this paper were carried out at LME in a wind-tunnel; atten- 
tion was mainly paid to the flow in the symmetry plane. The experimental device was adapted 
for low velocity flows with relatively small Reynolds numbers, that, nevertheless, provided tur- 
bulent regimes behind a backward-facing step. The results were obtained with a Particle Image 
Velocimetry (PIV) device. 

The numerical simulation of the flow was done on the basis of quasi-hydrodynamic (QHD) 
equations, proposed by Yu.V. Sheretov QHD equations differ from the Navier-Stokes system 
in additional dissipative terms of a second order in space. 



The scheme of the investigated flow is presented in Fig.^ The Reynolds number is given by 
Re = (UqK)/v, where Uq - is the average gas velocity in the entrance section of the channel, h - is 
the height of the step, v - is the kinematic viscosity of the medium. The first computational and 
experimental results were presented in d]-C7|. A brief description of the QHD system and its 
testing by computing a laminar backward- facing flow is given in the part I of the present publication 
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Figure 1: Wind tunnel section and notations. 



2 Experimental setup 

2.1 Wind tunnel 

The experimental results have been obtained in an open air-driven wind-tunnel at LME (Orleans, 
France). The existing square section (300 x 300mm) has been adapted to present a two-dimensional 
backward-facing step with adjustable step height. The width and the length (800 mm) of the test 
section were not modified. The electrically powered fan allowed us to obtain velocities of about 1 
to 80 m/s in the non-adapted test section without any obstacle. Previous works published in the 
literature about this subject |Tj pointed out that the reattachment length L of the flow behind the 
step varies with the Re number and the height h of the step: L can be twenty times larger than 
h. So, to be able to visualize the whole flowfield (and especially the reattachment area), the test 
section height H' upstream of the step was fixed to 51mm and we have considered different step 
heights between 12 mm and 50 mm. The length of the test section before the step is 200 mm, which 
allows the flow to be fully developed. The mean velocity obtained is approximately 1.4 m/s (as 
checked with a two-components hot-wire anemometer). The mean turbulent ratio then measured 
is less than 0.85% ahead of the step. One wall is transparent, for direct visualization, and the 
other ones are black. A laser sheet enters the wind tunnel section through a glass window in the 
upper wall. 

2.2 Particle Image Velocimetry (PIV) 

The experimental data generated are the two-dimensional components of the velocity behind the 
backward-facing step. An oil generator, placed at the entrance of the tranquillization room (velocity 
about 0.01 m/s), is used to generate and supply tracer particles. The mean diameter of the particles 
is about 1 fim. The laser sheet is generated by a double-oscillator laser: a Nd/Yag laser (Spectra 
Physics 400) adjusted on the second harmonic and emitting two pulses of 200 mJ each (A = 532 nm), 
at a repetition rate of 10 Hz. The laser sheet is developed with an optical arm containing mirrors. 
Lenses allow us to obtain a laser sheet with a divergence of about 60° and with a thickness of about 



1 mm in the vicinity of the step. For the present experimental data, the flow images are picked 
up by a P1VCAM CCD camera with 1008 x 1016 sensor elements, placed perpendicularly to the 
laser sheet. The laser pulses are synchronized with the image acquisition by a TSI synchronizer 
system driven by the InSight-NT™ software. As stated previously, the reattachment length can 
be as long as twenty times the height of the step. The aim of these experimental measurements 
is to obtain the mean flow behind the step, so we have chosen to decompose the flow area into 
several visualized areas (80 x 80 mm or 135 x 135 mm) to ensure a good precision in the wake of the 
step. This is only possible because the mean velocity field is studied. These sub-areas overlap by 
about 10 mm. For all the results presented here, the PIV recordings are divided into interrogation 
areas corresponding to 64 x 64 pixels. For data post-processing, the interrogation areas overlap 
by 50%. The local displacement vector is determined for each interrogation area by statistical 
methods (auto-correlation). The projection of the local flow velocity vector onto the laser sheet 
plane is calculated by InSight using the time delay between the two illuminations (At — 1 ms) 
and the magnification at imaging. The post-processing used here is very simple - no more than a 
velocity range filter. 

3 Numerical modeling of turbulent flows 

Numerical modeling of flow over a backward-facing step is carried out in accordance with the 
experimental parameters for the Reynolds numbers Re(h) = 4667, 4012, 1667 and the ratios 
h/H = 0.5, 0.44, 0.33 respectively. 
The flow in the wind-tunnel is air 

at room temperature, atmospheric pressure with entrance velocity Uq ~ 1.2 - 1.4 m/sec. The 
sound velocity in air under normal conditions equals c s = 340 m/sec. Because U$/c s ~ 0.003, we 
apply the approximation of a viscous incompressible isothermal flow |19|. Experimental visualiza- 
tion is done in the symmetry plane; therefore, we use a plane two-dimensional model of a flow in 
the numerical computations. 

These simplifying assumptions permit us to apply the finite difference algorithm described and 
tested in |18| and allowing for non-stationary flows. Unlike |18j . we now specify a plane velocity 
profile in the entrance section, that coincides with the conditions of the experiment. 

Numerical calculations were performed using the quasi-hydrodynamic (QHD) equations that 
differ from the Navier-Stokes ones by additional dissipative terms. These terms include a multi- 
plying factor r, - relaxation (smoothing) parameter, which has a dimension of time. For laminar 
flows, r is related to the molecular viscosity coefficient. For turbulent flows the value of r is no 
longer related with molecular viscosity and must be adjusted to fit the general flow features. The 
theory of QHD equations is presented in J3| and briefly described in [T8) . 

The results of computations are gathered in Tables ^CD an d are presented in Figs. [211761 
Additional computations may be found in |16| . Tables contain the serial number of computation 
(run), the dimcnsionless value of the smoothing parameter r, the spatial grid step, the number of 
grid nodes N y x N x , the dimcnsionless computing time To and the reference to the figures that 
illustrate this computation. Spatial grids are uniform and have equal steps in both directions 
(h x = h y ). Parameter r varies widely. The dimcnsionless time step in all presented variants is 
equal to St = 10 -4 . 

Velocity and pressure fields were recorded every 5000 steps, that is, with a time interval Ati = 
0.5. The velocity components in four particular locations behind the step were recorded every 500 
steps, i.e., with a time interval At — 500 St = 0.05. Energy spectra of the velocity components 
pulsations were computed for Re — 4667 and Re = 4012 according to the algorithm described in 
the Appendix. 

Computations were carried out in the dimcnsionless variables defined in |18|. For Re = 1667, we 
had H = 71 mm, Uo = 1.25 m/s, and the dimensionless time t = 1 corresponding to H/Uq = 0.056 s; 
for Re = 4012 and 4667, this variable corresponded 0.063 s and 0.072 s, respectively. 



3.1 Variant Re = 4667 

The major computational work for the proper choice of the smoothing parameter r was carried 
out for Re = 4667 (Tabd FigsHQ- 
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Figure 2: Instant stream functions (a), trajectories of fluid particles (b), averaged flow field (c) for 
Re=4667, run 2. 



Fig. demonstrates the sequences of the stream function isolines, computed according to 
algorithm at four successive times (indicated at the right of the figure. A single vortex 
appears behind the step; its length increases and then the single large vortex splits into several 
small ones that rotate in opposite directions. These smaller vortices are torn away from the main 
one and carried down the flow till they leave the computational domain. 

Fig. Eb demonstrates the trajectories of the fluid particles; the Fig. demonstrates the velocity 
field time-averaged between t\ and ti (values indicated in the figure). The averaged velocities 
and Uy" in every computational point are computed as 

< = T i T r u x (t)dt, = 

*2 - h Jt! t 2 -t 1 



/ Uy{t)dt. 
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Table 1: Computations for Re = 4667, H = 101mm, U = 1.40 m/s, h/H = 0.5. 
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Figure 3: Instant stream functions (a) and trajectories of fluid particles (b) for Re — 4667, run 4. 



The stream function is then computed using the averaged values. This procedure is correct due 
to the commutativity of temporal averaging and spatial differentiation. 

From Fig. [2] it is easily seen that for a quickly changing flow the flow pattern looks essentially 
different at four different times. The recirculation flow is well seen only on the stream functions 



that are constructed using the averaged velocities (Fig.^). The size of the separation zone equals 
L s /h ~ 3, which is twice less than observed experimentally. The same mean size of the separation 
zone is obtained in computations on a twice less accurate spatial grid with the same r = 0.001 (run 
3) and also with asmaller r = 0.0001 (run 1). In the last case the time step was equal to St = 10~ 5 . 
Therefore, for r < 0.001 a quasi-periodical flow is formed behind the step, and its shape rather 
weakly depends on the value of r and on the size of the spatial grid. In these computations the 
separation zone proves to be smaller than the one observed in experiments. 

Computations 4 and 5 correspond to r = 0.05. In Fig. [31 IH1 we demonstrate the results of 
computation 4. As well as in other computations, instantaneous pictures of the stream function 
(Fig. |3Ja) look much the same and strongly resemble those in Fig. The trajectories of fluid 
particles (Fig. |3Jd) do not exhibit a recirculation zone. 




a) 





Figure 4: Fragments of the averaged flow field: computational (a) and experimental (b). 

In FigQJi we show computed pictures of a flow fragment behind the step. This fragment 
corresponds to the region - 0.75 along the ordinate and 0.4 - 1.3 along the abscissa. Each numerical 
pattern is constructed from the instantaneous velocity vectors, averaged over time intervals t 2 —ti = 
0.015 (1 ms). We give 4 successive pictures starting from times T\ = 70, T 2 = 71, 5, T 3 = 73 and 
T4 = 74,5. So the numerical flowfields are obtained at a time interval of 1.5 (0.1 s, frequency 10 
Hz). 

FigEb shows the time sequence of experimental flowfields obtained at a frequency of 10 Hz in the 
same part of the flow domain behind the step. The velocity fields are obtained by intercorrelation 
between two images acquired with an interval time of 1 ms, that corresponds with the time delay 
between two laser illuminations. 

Averaging intervals in Figs0Ji and^Ja are approximately the same. Experimental and computa- 
tional portraits of the flow both demonstrate a chaotic non-stationary kind of flow, and the typical 
sizes of inhomogeneities obtained from computations agrees with the ones observed experimentally. 

Fig. [S] we show the computed flow fields, averaged over larger time intervals (£2 — t\ equals 
0.5 - 2 s). In Fig. [|)]we show the experimental flow fields averaged over intervals of 8-10 seconds. 
One can see the good agreement - both in size (L s /h ~ 6) and in structure of the separation 
zone in these averaged pictures. In both cases the averaged pictures depend on initial time t\ and 
integration time t 2 — t\, but this dependence is weak, because the averaging time is essentially 
larger than the typical time of velocity fluctuations. 

Fig. illustrates the time dependence of velocity components u x and u y as functions of time at 
three points located behind the step and far from the wall. The flow is essentially non-stationary 
and quasi-periodical. The general features of these curves depends weakly on the location of the 
points. 

The pulsation spectra E(k) of the velocity component u x at three spatial points are plotted in 
Fig. |SJ The mentioned spectra may be considered as reliable up to k k, 240. The k = harmonic 
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Figure 5: Averaged flow fields, computation, Re — 4667. 



is absent because it is related to the average velocity u. Fourier series expansion is carried out 
after completing the initial relaxation to stable oscillations during the time interval from t = 40 
to t = 120. The main oscillation frequencies and the decrease of amplitude for the harmonics at 
large and small values of k are easily seen. 

The results of computations 4 and 5 on different spatial grids were shown to be close to one 
another. 

Computation 6 is carried out with a large value of the smoothing parameter r. The structure 
of the flow proves to be non-physical - the non-stationary separation zone behind the step grows 
indefinitely and reaches the right border of the computation zone. In this case the velocities 
in the flow change slowly and the trajectories of the fluid particles are close to the streamlines. 
The averaged flow is of no interest in this computation. Additional computations and pictures 
illustrating them can be found in |16| . 




Figure 6: Averaged flow fields, experiment, Re = 4667. 




Figure 7: Fragment of temporal evolution of the Fi g ure 8: Velocity pulsation spectrum E{k), 
velocities, computation, Re = 4667. computation, Re = 4667. 

3.2 Variant Re = 4012 

Computation results for this variant in general repeat those obtained before. In run 1 (see Table 
2) we obtain a non-stationary solution with a separation zone length [L a /h ~ 2) smaller than the 
experimental one. 
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Table 2: Computations for Re = 4012, H = 92 mm, U = 1.45 m/s, h/H = 0.44. 



In computation 2 with t — 0.05 we get pictures of the flow averaged over ti—t\ = 0.015 (see 
Fig. In this figure, we show 8 successive pictures starting from T=I54 with a time interval 

1.5, that corresponds to a frequency of 10Hz. The corresponding experimental flow fragments 
are given in Fig. |5Jd starting from time t\ = 0.4 s with an acquisition rate equal also to 10Hz. 
Both starting points are chosen arbitrarily. Both calculated and experimental flow fields exhibit a 
chaotic non-stationary character. 




b) 

Figure 9: Fragments of the averaged flow field: computational (a) and experimental (b) for Re = 
4012. 

In Fig. E3 w e demonstrate the streamlines, averaged over long time intervals and constructed 
in the same way as in |18j . In Fig. 1111 the corresponding experimental pictures of the flow in the 
separation zone are shown over a background color that represents the velocity component u x . 

Fig. ^]depicts the time evolution of the velocity components behind the step. It demonstrates 
quasi-periodical oscillations. Figs. H3landH4lshow the energy spectrum of pulsations of the velocity 
components in linear and logarithmic scales, respectively. The last one contains the dependence 
E(k) ~ fc~ 5 / 3 (Kolmo gorov-Obukhov law in spectral form) This law of the kinetic energy 

dissipation is typical for well-developed turbulent flows. It can be easily seen that the numerical 
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Figure 10: Averaged flow fields, computation, Re = 4012. 
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Figure 11: Averaged flow fields, experiment, Re = 4012. 



energy spectrum in the frequency range, that agrees with the step of the spatial grid and with the 
time of computation To, coincides in general with the classical law of the decreasing of energy E 
with the growth of the wave number k. 
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Figure 12: Temporal evolution of the velocity behind the step, computation, Re = 4012. 
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Figure 13: Velocity pulsation spectrum E(k), Fi g ure 14: Velocity pulsation spectrum E(k), 
computation Re = 4012 Logarithmic scale, computation, Re = 4012. 



3.3 Variant Re = 1667 

For Re — 1667 (see Table 3), in run 2, the length and the structure of the separation zone 
(L s /h ~ 6) agrees with the experimental data. In run 1 (L s /h ~ 2, 5) the length of the separation 
zone proves to be smaller than in the experiment. In run 3 the numerical solution becomes non- 
physical - the separation zone grows indefinitely. 

In Fig. 115b we show a series of the averaged stream function fields (run 2) The corresponding 
averaged experimental flow picture is shown in Fig. 115b . The size and the structure of the flow 
in the separation zone are in good agreement. Fig. 1161 illustrates the temporal evolution of flow 
velocity components; it presents a quasi-periodical character. 
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Table 3: Computations for Re = 1667, H = 71mm, U = 1.25 m/s, h/H = 0.33. 




Figure 15: Averaged flow fields: computational (a) and experimental (b) for Re = 1667. 



4 Discussion and conclusions 

4.1 Relaxation parameter 

Numerical computations are based on the QHD equations. They differ from the Navier-Stokes 
system by additional dissipative terms that contain a small parameter r - relaxation (smoothing) 



Re = 1667 




100 110 120 130 140 150 

Time 



Figure 16: Temporal velocity evolution behind the step, computation, Re = 1667. 



parameter. While computing turbulent flows, r is chosen much larger than for computing laminar 
flows |18) . We manage to obtain the structure of the flow behind a backward-facing step by varying 
r. The dependence of the averaged flow on t proves to be relatively weak: in the variant Re — 4667 
the length of the separation zone varied from L s /h ~ 3 to L s /h ~ 6, when r was varied from 0.001 
to 0.05. We should also mention that the parameter of relaxation r is the only constant in the QHD 
model that must be adjusted to achieve agreement between numerical and experimental results. 

A posteriori specification of the relaxation parameter on the basis of the flow field resembles 
to adjusting the turbulent viscosity in phenomenological models of turbulence. In these models 
the value of turbulent viscosity is not related with the coefficient of the molecular viscosity. This 
quantity is derived from experimental data, e.g., from the pulsation features of the flow or from 
its characteristic scales |19) . 

Agreement with the experimental data - as far as the length and the structure of the separation 
zone are concerned - is obtained in variants Re — 4667 and Re — 4012 for r = 0.05, and in 
variant Re — 1667 for t = 0.02. In dimensional form it corresponds to 3.6 • 10~ 3 s, 3.2 • 10~ 3 s and 
1.12 • 10 -3 s respectively. In these cases the average length of the separation zone equals L s /h ~ 6, 
and the position of the vortex centers coincide in computations and in experiments. The position 
of a secondary vortex located at the internal corner is also well described in the computations. 

Mention, that in our calculations we have not seen the well-known anomaly of low-Reynolds- 
number (LRN) Reynolds-stress models (RSM) (|7|, [S]), that is a spurious bending (backward 
curvature) of the dividing streamline at the reattachment point. 

There remains the problem of a priori choice of the relaxation parameter in QHD equations for 
computing turbulent flows. Nevertheless, we may claim that r < 1, because the dissipative terms 
should be considerably smaller than the convective ones. In computations the minimal resolution 
scale is the step of the spatial grid. Therefore we may suppose that r should be greater than the 
characteristic perturbation spreading time over one spatial step h x ; that is, r > h x /c s . Besides, 
for the problem under consideration, r should be a fortiori greater than the averaging time, for 
which the velocity profile at the entrance of the channel may be considered as a flat. This time 
was not measured in the present experiments. 



4.2 Flow visualization and time averaging intervals 

The numerically observed type of the flow formation resembles the well-known Landau turbulence 
scenario ^5]. This formation consists of the growth of initial perturbations that are suppressed by 
nonlinear media reactions. It leads to the appearance of a complex quasi-periodical motion, which 
seems to be turbulent [JCL6]. 

The computational time step St that describes the variation of the nonstationary velocity 
field, is several orders of magnitude smaller than the time resolution of the experiments. So it is 
impossible to compare the numerical and experimental flow patterns directly, using instantaneous 
streamlines. The trajectories of the liquid particles, based on the stream- functions, are not suitable 
for the comparison with the present experiments. That is why we used averaged velocity fields 
when comparing numerical and experimental results. 

We used two types of temporal intervals for visualizing the results of computations - large 
intervals for comparing global flow pictures and small intervals for comparing local " snapshots" . 

The averaging interval for global pictures equals t 2 — t\ = 10 — 50 dimensionless units, that is, 
0.5-3 s. This time interval exceeds largely the averaged period of the temporal velocity evolution 
(Figs. 171 IT21 and lPo|) . The number of averaged fields proves to be large enough, because the flow 
recording is done with the interval 0.5. In other words, the number of averaged velocity fields in 
computation approximately equals the number of averaged experimental images. It makes both 
averaging procedures equivalent in general. The obtained fields depend weakly on the beginning 
and length of the time averaging interval. 

The averaging interval for comparing numerical and experimental "snapshots" equals to 0.015 
dimensionless units, approximately 1 ms, which coincides with the time delay between the two 
illuminations - a gap between two pulses of the PIV installation. The numerical flowfields are ob- 
tained at a time interval ~ 0.1 s., to reproduce the experimental acquisition frequency of frequency 
~ 10 Hz. The computational and the experimental averaged pictures in general look alike - both 
consist of small chaotically mixed non-stationary whirls that essentially depend on the choice of 
the temporal interval t 2 —t\. 

It is impossible to observe instant velocity fields experimentally. Therefore the question of 
validity of the computed instant velocity fields is physically senseless. At the same time numerical 
and experimental fields averaged over large and small time intervals prove to be consistent with 
each other. 

4.3 Conclusions 

In |18j the QHD equations were tested in a laminar flow simulation. There a relaxation parameter 
r was determined by the molecular viscosity. In computations it was chosen as small as possible 
to provide the stability of the numerical computations. The stationary flow did not depend on its 
value. The initial oscillations of the solution faded during the time. 

In the present paper we demonstrate that in the turbulent flow the initial perturbations do not 
fade and lead to the formation of non-stationary flow pattern, that resembles a periodical process. 

This type of solution is observed for sufficiently small values of the dimensionless temporal 
relaxation parameter: r < 0.1. 

The comparison of numerical solutions with experimental data is possible only on the basis of 
the averaged velocity fields. In this case the outlook of the averaged flow depends both on the value 
of the temporal relaxation parameter r (which makes the difference between the QHD system and 
the Navier-Stokes equations) and on the chosen averaging interval. 

In spite of considerable simplifications of the problem - using a two-dimensional model of the 
viscous non-compressible isothermal fluid - the results of computation are in good agreement with 
the present experiments. 

We suppose that the numerical modeling of the turbulent flow over a backward-facing step by 
means of the QHD system is possible due to the additional (in comparison with the Navier-Stokes 
equations) dissipative terms that are small for the stationary flows and become not small for the 
non-stationary turbulent flows. 



5 Appendix. Computing the spectrum of pulsations 

The energy spectrum of the flow pulsations is computed on the basis of the time dependence of 
velocity components u x and u y . Let us write the velocity u x as 



where u x is the average velocity value, u' x is its pulsation component. Then we apply the discrete 
Fourier transform to u' x (for simplicity further on we omit primes). We present the pulsation 
component in form of a series 



where u x - is the value of the velocity u x at i,j at time I At; N is the number of the considered 
records, k runs from 1 to N/2. 

The energy spectrum is derived from these Fourier coefficients: 



The maximal value of k = N/2 corresponds to oscillations of period At. When the frequency of 
recordings is chosen, oscillations with higher frequencies are not allowed. But for proper resolution 
we need at least 4-5 points per period. Therefore the obtained spectrum can be considered as 
reliable for k < N/8. 

For all variants of the present calculations, At — 0.05. For example, for Re = 4012, run 2, 
Tq = 200, so the reliable range of k in Figs. fTSlHH is limited to k < T /At/8 = 500 approximately. 
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